As pointed out in the community about the hidden cost of median(x) and possible fix.

Summary

This change replaces the full sort used in median() for dense inputs with order-statistic selection via nth_element, reducing average complexity from O(n log n) to O(n). Sparse inputs retain the original sort-based code path to preserve historical behavior and output sparsity.

Algorithmic change

Previous implementation:

• sort(x) followed by indexing middle element(s)

New implementation:

• nth_element(x, k) for odd length
• nth_element(x, [k, k+1]) for even length

This computes only the required order statistics without fully sorting.

NaN handling

NaN semantics are preserved explicitly:

• "omitnan": NaNs are filtered before selection.
• "includenan": slices containing NaNs return NaN.

This matches the behavior previously provided implicitly by sort()

Sparse behavior

Sparse inputs continue to use the original sort-based code path to preserve:

• sparsity of outputs,
• historical indexing behavior,
• and memory characteristics.

Dense inputs use the new nth_element path.

I went a step ahead and ran a benchmark to see if it's really doing the speed up as promised, here are the results below. I am also attaching the benchmarking script.

Performance

Benchmarks on large dense inputs (Octave 9.x, Apple clang):

Speedup grows with input size, consistent with removal of the log n factor.

Benchmarking Script can be found here
CSV results from which plots were generated can be found here

Notes

• No semantic or API changes.
• Only the internal selection algorithm is modified.
• All existing tests pass
• All existing BISTs pass unchanged.
• As mentioned in community, NaN values are taken care of properly
• I can move to improving the implementation of KDTree and improve it's build time if this gets merged as median was being a bottleneck as raised in this issue